{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "import scipy.linalg as la\n",
    "import scipy\n",
    "from scipy.linalg import expm, sinm, cosm\n",
    "\n",
    "import sympy\n",
    "from sympy import Matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = sympy.Symbol(\"h\",real= True)\n",
    "P = Matrix([ [-1,1,0],\n",
    "               [0,-1,1],\n",
    "               [0,0,1]   ])\n",
    "D = Matrix([ [-2,1,0],\n",
    "               [1,-2,1],\n",
    "               [0,1,-2]  ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# matrix P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = \n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 1 & 0\\\\0 & -1 & 1\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1,  1, 0],\n",
       "[ 0, -1, 1],\n",
       "[ 0,  0, 1]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"P = \")\n",
    "P\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.is_diagonal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.is_diagonalizable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1, Ja = P.jordan_form()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 1 & 0\\\\0 & -1 & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1,  1, 0],\n",
       "[ 0, -1, 0],\n",
       "[ 0,  0, 1]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & \\frac{1}{4}\\\\0 & 1 & \\frac{1}{2}\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1, 0, 1/4],\n",
       "[0, 1, 1/2],\n",
       "[0, 0,   1]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 1 & 0\\\\0 & -1 & 1\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1,  1, 0],\n",
       "[ 0, -1, 1],\n",
       "[ 0,  0, 1]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P1 * Ja * P1.inv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{h}$"
      ],
      "text/plain": [
       "exp(h)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = sympy.E**h\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- e^{h} & e^{h} & 0\\\\0 & - e^{h} & e^{h}\\\\0 & 0 & e^{h}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-exp(h),  exp(h),      0],\n",
       "[      0, -exp(h), exp(h)],\n",
       "[      0,       0, exp(h)]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = P1 * (Ja*f) * P1.inv()\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matrix D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 1 & 0\\\\1 & -2 & 1\\\\0 & 1 & -2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-2,  1,  0],\n",
       "[ 1, -2,  1],\n",
       "[ 0,  1, -2]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D.is_diagonalizable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q,D1=D.diagonalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 0 & 0\\\\0 & -2 - \\sqrt{2} & 0\\\\0 & 0 & -2 + \\sqrt{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-2,            0,            0],\n",
       "[ 0, -2 - sqrt(2),            0],\n",
       "[ 0,            0, -2 + sqrt(2)]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 1 & 1\\\\0 & - \\sqrt{2} & \\sqrt{2}\\\\1 & 1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1,        1,       1],\n",
       "[ 0, -sqrt(2), sqrt(2)],\n",
       "[ 1,        1,       1]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{1}{2} & 0 & \\frac{1}{2}\\\\\\frac{1}{4} & - \\frac{\\sqrt{2}}{4} & \\frac{1}{4}\\\\\\frac{1}{4} & \\frac{\\sqrt{2}}{4} & \\frac{1}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1/2,          0, 1/2],\n",
       "[ 1/4, -sqrt(2)/4, 1/4],\n",
       "[ 1/4,  sqrt(2)/4, 1/4]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q.inv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 2 e^{h} & 0 & 0\\\\0 & \\left(-2 - \\sqrt{2}\\right) e^{h} & 0\\\\0 & 0 & \\left(-2 + \\sqrt{2}\\right) e^{h}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-2*exp(h),                     0,                     0],\n",
       "[        0, (-2 - sqrt(2))*exp(h),                     0],\n",
       "[        0,                     0, (-2 + sqrt(2))*exp(h)]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D1 * sympy.E**h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- e^{h} + \\frac{\\left(-2 - \\sqrt{2}\\right) e^{h}}{4} + \\frac{\\left(-2 + \\sqrt{2}\\right) e^{h}}{4} & \\frac{\\sqrt{2} \\left(-2 + \\sqrt{2}\\right) e^{h}}{4} - \\frac{\\sqrt{2} \\left(-2 - \\sqrt{2}\\right) e^{h}}{4} & \\frac{\\left(-2 - \\sqrt{2}\\right) e^{h}}{4} + \\frac{\\left(-2 + \\sqrt{2}\\right) e^{h}}{4} + e^{h}\\\\\\frac{\\sqrt{2} \\left(-2 + \\sqrt{2}\\right) e^{h}}{4} - \\frac{\\sqrt{2} \\left(-2 - \\sqrt{2}\\right) e^{h}}{4} & \\frac{\\left(-2 - \\sqrt{2}\\right) e^{h}}{2} + \\frac{\\left(-2 + \\sqrt{2}\\right) e^{h}}{2} & \\frac{\\sqrt{2} \\left(-2 + \\sqrt{2}\\right) e^{h}}{4} - \\frac{\\sqrt{2} \\left(-2 - \\sqrt{2}\\right) e^{h}}{4}\\\\\\frac{\\left(-2 - \\sqrt{2}\\right) e^{h}}{4} + \\frac{\\left(-2 + \\sqrt{2}\\right) e^{h}}{4} + e^{h} & \\frac{\\sqrt{2} \\left(-2 + \\sqrt{2}\\right) e^{h}}{4} - \\frac{\\sqrt{2} \\left(-2 - \\sqrt{2}\\right) e^{h}}{4} & - e^{h} + \\frac{\\left(-2 - \\sqrt{2}\\right) e^{h}}{4} + \\frac{\\left(-2 + \\sqrt{2}\\right) e^{h}}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[      -exp(h) + (-2 - sqrt(2))*exp(h)/4 + (-2 + sqrt(2))*exp(h)/4, sqrt(2)*(-2 + sqrt(2))*exp(h)/4 - sqrt(2)*(-2 - sqrt(2))*exp(h)/4,        (-2 - sqrt(2))*exp(h)/4 + (-2 + sqrt(2))*exp(h)/4 + exp(h)],\n",
       "[sqrt(2)*(-2 + sqrt(2))*exp(h)/4 - sqrt(2)*(-2 - sqrt(2))*exp(h)/4,                 (-2 - sqrt(2))*exp(h)/2 + (-2 + sqrt(2))*exp(h)/2, sqrt(2)*(-2 + sqrt(2))*exp(h)/4 - sqrt(2)*(-2 - sqrt(2))*exp(h)/4],\n",
       "[       (-2 - sqrt(2))*exp(h)/4 + (-2 + sqrt(2))*exp(h)/4 + exp(h), sqrt(2)*(-2 + sqrt(2))*exp(h)/4 - sqrt(2)*(-2 - sqrt(2))*exp(h)/4,       -exp(h) + (-2 - sqrt(2))*exp(h)/4 + (-2 + sqrt(2))*exp(h)/4]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q * (D1 * sympy.E**h) * Q.inv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High Dimention"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0.],\n",
       "       [0., 0., 0.]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix([ [1,2,3],[4,5,6]]) #,[7,8,9]\n",
    "m,n=A.shape\n",
    "np.zeros(A.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "x = sympy.Symbol(\"x\")\n",
    "f = sympy.Function(\"f\")(x)\n",
    "def Jordan_block_F(Ja):\n",
    "    B = np.zeros(A.shape)\n",
    "    m ,n = A.shape\n",
    "    for i in range(m):\n",
    "        a = 0\n",
    "        for j in range(i,n):\n",
    "            g=f.diff(x,a)                                    #这里是错的 应该是Jordan块的阶次  先分解为约旦块才行\n",
    "            x = A[i,j]\n",
    "            B[i,j] = g(x)\n",
    "            a+=1\n",
    "    B = Matrix(B)\n",
    "    return B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# test others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# expm( np.array([ [2,0],[0,3] ]) )\n",
    "# expm( np.array([ [2,0],[0,3] ] * h ))  \n",
    "# expm( np.array( [ [2,0],[0,3]  ])* h )  \n",
    "# expm( h * Matrix( [ [2,0],[0,3]  ])  ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'expm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-22-16c67f1ee343>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# expm( Matrix([[2,0],[0,3]]) )  # wrong\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;31m# expm( h* np.array([[2,0],[0,3]]) )\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mexpm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'expm' is not defined"
     ]
    }
   ],
   "source": [
    "# expm( Matrix([[2,0],[0,3]]) )  # wrong\n",
    "# expm( h* np.array([[2,0],[0,3]]) ) \n",
    "expm( np.array([[2,0],[0,3]]) ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'expm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-23-1136263c546a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mexpm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mMatrix\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m  \u001b[1;31m# wrong, must array\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'expm' is not defined"
     ]
    }
   ],
   "source": [
    "expm( Matrix([[2,0],[0,3]]) )  # wrong, must array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'expm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-24-ebfc87c11f2c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msympy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSymbol\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"h\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mexpm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mh\u001b[0m\u001b[1;33m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'expm' is not defined"
     ]
    }
   ],
   "source": [
    "h = sympy.Symbol(\"h\")\n",
    "expm( h* np.array([[2,0],[0,3]]) ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'expm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-25-33290c0838da>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m3\u001b[0m \u001b[1;31m# h is symbol\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mexpm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mh\u001b[0m\u001b[1;33m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m  \u001b[1;31m# can be solve numerically\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'expm' is not defined"
     ]
    }
   ],
   "source": [
    "h = 3 # h is symbol\n",
    "expm( h* np.array([[2,0],[0,3]]) )  # can be solve numerically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[6, 0],\n",
       "       [0, 9]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h* np.array([[2,0],[0,3]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}6 & 0\\\\0 & 9\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[6, 0],\n",
       "[0, 9]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix([[2,0],[0,3]])\n",
    "h*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unsupported operand type(s) for ** or pow(): 'Exp1' and 'MutableDenseMatrix'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-28-905535eab1dc>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msympy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mE\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mh\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: unsupported operand type(s) for ** or pow(): 'Exp1' and 'MutableDenseMatrix'"
     ]
    }
   ],
   "source": [
    "f = sympy.E**(h*A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[exp(6), 1],\n",
       "       [1, exp(9)]], dtype=object)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[2,0],[0,3]])\n",
    "f = sympy.E**(h*A)\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# scipy idea"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as la\n",
    "import scipy\n",
    "from scipy.linalg import expm, sinm, cosm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np \n",
    "from scipy.linalg import expm, sinm, cosm\n",
    "expm(np.zeros((2,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Euler’s identity $ (exp(i*\\theta) = cos(\\theta) + i*sin(\\theta)) $ applied  to a matrix: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.4264593 +1.89217551j, -2.13721484-0.97811252j],\n",
       "       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([[1.0, 2.0], [-1.0, 3.0]])\n",
    "expm(1j*a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.4264593 +1.89217551j, -2.13721484-0.97811252j],\n",
       "       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cosm(a) + 1j*sinm(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 自由 矩阵函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4., 15.],\n",
       "       [ 5., 19.]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.linalg import funm\n",
    "import numpy as np \n",
    "a = np.array([[1.0, 3.0], [1.0, 4.0]])\n",
    "funm(a, lambda x: x*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4., 15.],\n",
       "       [ 5., 19.]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a.dot(a)\n",
    "a @ a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import eigh\n",
    "def funm_herm(a, func, check_finite=False):\n",
    "    w, v = eigh(a, check_finite=check_finite)\n",
    "    ## if you further know that your matrix is positive semidefinite,\n",
    "    ## you can optionally guard against precision errors by doing\n",
    "    # w = np.maximum(w, 0)\n",
    "    w = func(w)\n",
    "    return (v * w).dot(v.conj().T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.3890561 ,  0.        ],\n",
       "       [ 0.        , 20.08553692]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expm( np.array([[2,0],[0,3]]) ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2371846.83476555, 1081759.23958229],\n",
       "       [7572314.67707605, 3453606.07434784]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = 3\n",
    "expm( x * np.array([[2,1],[7,3]]) )   # Can be solved numerically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "2a6cc8b4fade00e9cef0cb90e5519d8e8bd0936e8153a2638b9472c580ba270d"
  },
  "kernelspec": {
   "display_name": "Python 3.8.8 64-bit ('base': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
